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Abstract 

The Bose-Hubbard model exhibits a rich phase diagram consisting both of 
insulating regimes where diagonal long range (solid) order dominates as well 
as conducting regimes where off diagonal long range order (superfluidity) is 



'^ ' present. In this paper we describe the results of Quantum Monte Carlo cal- 



culations of the phase diagram, both for the hard and soft core cases, with a 
particular focus on the possibility of simultaneous superfluid and solid order. 



H ' We also discuss the appearance of phase separation in the model. The sim- 

ulations are compared with analytic calculations of the phase diagram and 
spin wave dispersion. 
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I. INTRODUCTION 

A lot of attention has been focussed on the interacting electron problem in the last 
several decades, whereas the interacting boson problem has been considered more often in 
the framework of specific applications only. However, there are a number of important 
situations where the elementary excitations are either intrinsically bosonic in character or 
else can usefully be viewed in terms of bosonic models. '^He is an example of the former 
situation, |T]] while quantum spin systems, ^ granular superconductors, |^ and flux lines 
in type-II superconductors |^ are examples of the latter. Therefore it is important to 
understand in detail the features of model boson systems, in much the same way that one 
studies the Hubbard, Anderson, and t-J Hamiltonians for correlated fermions. In this paper 
we consider a lattice model of interacting bosons, the Bose Hubbard (BH) Hamiltonian: 

H = -t ^(ajaj + a]ai) - fij^^i + Vo^n^ + ViY^ riiUj + V2 X! ^i^k ■ (1) 

{ij) i i (ij) {{ik)) 

Here a, is a boson annihilation operator at site i, and rij = ajai. The transfer integral t = 1 
sets the scale of the energy, and /i is the chemical potential. Vo,Vi, and V2 are on-site, 
near-neighbor, and next-near-neighbor boson-boson repulsions. 

The interactions Vq, Vi, and V2 promote the formation of "solid" order, where the boson 
occupations fall into regular patterns, at special densities commensurate with the lattice. 
The hopping matrix element t favors mobile bosons, and consequently a superfluid phase at 
T = 0. In what follows the nature of the correlation functions will be studied as we change 
the Hamiltonian parameters and the density p = j^ Z]i(^i)- 

When Vq = cxD, the BH model maps onto the quantum spin-1/2 Hamiltonian 

H = -t Y^istsj + spn + v^y: sts^ + 1^2 E St SI -H,J2 st. (2) 

{tj) {ij) {(ik)) i 

The field H^ = p — 2Vi — 2V2. Since rii <-* S^ + |, ordering of the density corresponds 
to finite wave vector Ising type order. Similarly, a, ^^ S^^ so that superfiuidity maps to 
ferromagnetic ordering in the XY plane. One of the things we shall be interested in in this 



work is the possibility that density and superfluid order are not mutually exclusive. Indeed, 
aX V2 = Hz = 0, the special point Vi = 2t corresponds to the Heisenberg Hamiltonian where 
Ising and XY order coexist. It has been suggested by various authors [|,0 that the addition 
of further terms like V2 or Hz could stabilize this "supersolid" from a special symmetry point 
to a broader area of the phase diagram. Precisely at the Heisenberg antiferromagnet (AF), 
the effect of a field H^ is known: It breaks the full rotational symmetry and selects ordering 
in the XY plane since the spins can more easily take advantage of the field energy. This 
argument has been used to suggest why doping favors the superconducting over the CDW 
state in the negative-f/ Hubbard model where an analogous "supersolid" symmetry exists 
at half-filling. 

While there have been many mean field (MF) studies of the spin Hamiltonian Eq. 2, there 
have been to date only a few numerical studies PHTOH- Monte Carlo studies of interacting 



quantum boson and spin models provide a useful, exact method to study the nature of 
the correlations on finite lattices. Combined with finite size scaling methods, they can be 
used to extract information concerning the thermodynamic limit. Boson simulations are 
somewhat easier than related path integral methods for interacting electron systems, since 
they can utilize algorithms which scale linearly with the lattice size and can reach essentially 
arbitrarily low temperatures. 

This paper is organized as follows: In Section II and HI we determine analytically the MF 
phase diagram, extending past work by considering additional types of order, and describe 
spin wave calculations of the dispersion relations in the various phases. In Section IV 



we provide numerical results for the soft core model, extending our earlier studies [jTO 
In Section V we describe new results for the hard-core phase diagram. Conclusions are 
presented in Section VI. 



II. MEAN FIELD PHASE DIAGRAM 



Previous work established the MF phase diagram of the spin Hamiltonian considering 
only the possibility of superfluidity and Neel-type ordering of the density. |^,^.11| At half- 



filling, or equivalently at zero magnetization M^ = 0, for V\ > 2 and < V2 < Vi — 2 
the spins form a Neel state, corresponding to a checkerboard Bose solid with an ordering 
vector k* = (7r,7r). For V2 > max (Vi — 2; 0) a ferromagnetic phase is formed, with a net 
moment M^y 7^ and M^ 7^ 0. This phase corresponds to a superfluid, and is also stable 
for arbitrary Vi and V2 away from half-filling. A fully polarized magnetic phase in a strong 
magnetic field Hz, where only Mz 7^ 0, corresponds to a Mott-insulator with precisely one 
boson per site. As the solid and the superfluid phases possess different broken symmetries, 
one could expect that the transition between them is first order. However, a rather different 
scenario has also been put forward, suggesting that the - presumably - first order transition 
is split up into two distinct second order transitions, where the two order parameters vanish 



at separate points |[TT|JT^jr5[| In the regime between the two transitions both order parameters 
are non-zero, hence it has been termed a supersolid. |]5|,|6|,pJ^ This intriguing possibility is 
the subject of the investigations reported in this paper. 



The mean field analysis indeed finds such a supersolid phase P,PJTl[], although in the 
hard core limit longer range forces (V2 > 0) are needed to stabilize it. However, recently 
it was claimed that this conclusion changes in the soft core case, and a supersolid phase 



exists with nearest neighbor interaction alone [1^. Finally, recent studies on the related 
Heisenberg model with competing first and second neighbor couplings Ji and J2 established 
the possibility of additional phases: a collinear phase, with alternating lines of up and down 



spins, at large J2/J1 [OlHl' ^^^ ^ disordered phase at intermediate values of J2/ Ji- pO|j2T 
These differing results clearly call for a reinvestigation of the problem. 

The MF phase diagram of the spin Hamiltonian Eq. 2 worked out by Matsuda and 
Tsuneto 0, and described above, allowed only for a two-sublattice magnetic ordering of 
the spins corresponding to a Neel solid. Representing the spins by classical vectors of length 



S we extend earlier MF analyses pJT^ for the case of a square lattice to include also the 
possibility of a coUinear phase which is expected to form for intermediate to large next- 
near-neighbor repulsion V2 (see Fig. 1). Assuming that the spins are ordered in the XZ 
plane, the MF energies per spin, cn and ec, of the Neel and collinear spin configurations 
are given by 

e^ = -45*^ sin 6 a sin Ob + 25" Vi cos 6 a cos Ob + V^25'^ (cos^ 6a + cos^ Ob) - 

TT 

— ^Sycos9A + cosObj (3) 

T/ Q2 

ec = -S (sin6'/ji + sin6'ij2) H — {cos9ri + cos9r2) + 2V2S cos9ri cos9r2 - 

IT 

— ^S {cos9ri + cos9r2) . (4) 

9a and 9b are the angles between the spin direction and the 2;-axis on sublattice A and B, 
respectively. 9ri and 9r2 are the corresponding angles in the collinear phase on even and 
odd rows. The different phases are identified as follows: 

cos 9 A = cos 9b < 1 or cos9ri = cos 9r2 < 1 Superfluid 

cos 9a = — cos 9b = I Neel Solid 

cos 9ri = — cos 9r2 = 1 Collinear Solid 

sin 9a 7^ sin 9b and — 1 < cos 9a 7^ — cos 9b < I Neel Supersolid 

sin 9ri 7^ sin 9r2 and — 1 < cos 9ri 7^ — cos 9r2 < 1 Collinear Supersolid 

cos 9a = cos 9b = i or cos 9ri = cos 9r2 = 1 Mott phase (5) 

We performed the MF analysis in the same spirit as in Refs. ||5|,|6|JTl[l . One proceeds by 
minimizing cn and ec separately with respect to the angles 9a, 9b and 9ri, 9r2, respectively. 
Then the results for fixed magnetic field Hz are translated to fixed magnetisation, i.e. boson 
density. Finally, we compare the energies of the different phases to obtain the complete 
MF phase diagram of the spin Hamiltonian Eq. 2. Explicitly, for two-sublattice Neel type 
ordering we find the following phases for Vi > 2 and < V2 < V^i — 2: 

m= Solid 



1 /Vl - ^2 - 2 

< m< -\ Yj 77 Neel Supersolid 

2 V Vi — V2 + 2 



1 l^l-\/2-2 1 CO., 

2\l v,-V, + 2 ^'^^2 Superfluid 

171= - Mott Insulator (6) 

where m = |p— || is the magnetisation of the system. For < Vi < 2 there is no Neel order 
and for m ^ the MF ground state is always superfluid. 

Similarly we analyze the phase diagram following from minimizing ec for the ordered 
collinear spin structures corresponding to an ordering wave vector k* = (0, tt) or (vr, 0). 
At half-filling the collinear solid (see Fig. 1) is realized for arbitrary values of the near- 
neighbor repulsion Vi. The reason is that at half-filling the energy for the collinear solid 
is ec = — V2/2, i.e. independent of Vi due to the cancellation of S^S^ energies for near- 
neighbor sites on the same and neighboring rows. Away from half-filling only the superfluid 
minimizes ec for V2 < 2. For V2 > 2 a collinear supersolid appears in the phase diagram 
and the boundary between the superfluid and the collinear supersolid is determined by 



1 /V'2-2 
< m< -\ Collinear Supersohd 

2V V2 



M^2-2 1 on., 

-W— ^^<m<- Superfluid (7) 

which is again independent of Vi. For V2 > 2 the collinear supersolid phase occurs in a 



density strip of width J{V2 — 2)/V2 around half-filling. 

Given the MF solution for e^ and ec separately, a comparison for the energies of the 
different phases allows to map out the complete mean-field phase diagram of the spin Hamil- 
tonian Eq. 2. E.g. at half-filling, m = 0, we have to compare 

ec = — V2 Collinear Solid 

ctv = ec7 = — 1 Superfluid 

ctv = -{V2- Vi) Neel Sohd (8) 

The resulting phase diagram is shown in Fig. 2. Interestingly, for 2 < Vi < 4 increasing V2 
drives two transitions: first increasing V2 frustrates the Neel solid and leads to a transition 

6 



to a superfluid. Increasing V2 further stabilizes collinear order and leads to a transition from 
a superfluid to a collinear solid. pT| 

Away from half-filling, < m < 1/2, no solids, neither Neel nor collinear are MF 
solutions. Instead, transitions occur between the superfluid, and the Neel and collinear 
supersolid phases. The boundaries between the different phases are given by 

V2 = Vl — 2 — - Superfluid to Neel Supersolid 

1 — 4m^ 

2 
V2 = 7^ Superfluid to Collinear Supersolid 

1 Am? 

V2 = -V^i — 7^ Neel to Collinear Supersolid (9) 

2 1 — 4m^ 

Finite doping leads to a rigid shift of the phase boundary lines obtained at half-filling with 
the solid replaced by supersolid phases. For 

2i±±!^ < V, < 4i±^ (10) 

1 - Am? 1 - 4m2 ^ ^ 

this still allows for two transitions with increasing V2, from a Neel supersolid to a superfluid 
to a collinear supersolid. The V\-V2 phase diagram for a fixed magnetisation m = 0.2 is 
shown in Fig. 3. In addition. Figs. 4 and 5 show the phase boundaries in the V2, mi plane 
for a fixed value of Vi and in the Vi, m plane for a fixed value of V2, respectively. 

Recently it was claimed that a finite core repulsion Vq < 00 qualitatively changes this 
picture. |T^ Supersolids were found to exist even at half-filling, moreover without the next 
nearest neighbor repulsion V2. To study these claims we extend the MF analysis by introduc- 
ing an approximate soft core representation allowing the spin length 5 to be a variational 
parameter and adding a term Hconstraint = VoZ]i(S^ — 1)^ to the Hamiltonian. The min- 
imization of the ground state energy is now done separately with respect to S"^, S^ and 

cA cB 

We expand the ground state energy around the superfluid phase, and consider the eigen- 
values corresponding to small spatial modulations of the density and superfluid order pa- 
rameter, in effect generating a Ginzburg-Landau type expression. The superfluid-collinear 
supersolid transition is studied by writing 



S^ = m — e S^ = m + e 



^ z 



S^ = s-6 S^ = s + 6 (11) 

and expanding to second order in the (small) fluctuations e and 6. The expectation value e 
of the ground state energy per site takes the form 

e = esF- 41/2e^ + Vq [(125^ + Am^ - 1)5^ + {I2m^ + As"^ - l)e^ + I6sme5] , 
esF = -8s2 + 4(\/i + V2)m^ + lvo{4s^ + 4m' - 1)' . (12) 

o 

The ground state energy is the sum of eigenvalues of a matrix in the (e, 6) space. First we 
solve for s at fixed number of particles, i.e. fixed m, in the superconducting state where 
S = e = 0, and obtain 

.^ = i-m= + A (13) 

A zero eigenvalue of the energy matrix signals the phase transition. The condition for the 
vanishing of the determinant can be solved for V2 for arbitrary m 

which gives the phase boundary between the superfiuid and the coUinear supersolid. With 
the same procedure the phase boundary between the superfiuid and the Neel supersolid is at 
V2 = Vl — 2 — 16m'/ [1 — Am"^ + 16/ Vq]. As in the hard core case the phase diagram displays 
Neel- and coUinear supersolid, and superfiuid phases. At half-filling the supersolid phases 
vanish, and two insulating solids are direct neighbors to the superfiuid, in contrast with the 
result of Ref. 0. This result is independent of Vq, i.e. it is true both in the soft and hard 
core limits, in agreement with the above hard core MF calculation. 

III. SPIN WAVE ANALYSIS 

The analyses of the spin wave fluctuations which exist in the literature P, pT]JTB| are in 
disagreement. The spectrum has been found to be either linear |^ or quadratic |lTTl , p!3| at 



the solid-supersolid phase boundary. This dependence is crucial for numerical studies, as it 
determines the dynamical critical exponent z and thereby the appropriate finite size scaling 
of the lattice. 

To settle the issue, we redo the linear spin wave theory analysis for the spin model of 
Eq. 2 and determine the spectrum in the superfiuid, the Neel solid and the Neel supersolid. 
Again we assume that the spins are ordered in the XZ plane with an angle Qa{B) to the 
2;-direction. On each sublattice the spin quantisation axis is rotated to align the spins along 
the local direction of the magnetisation by 



'iGA(ie-B) 



COS Qa{b) - sin Qa{b) 



^ieA(jeB) 



(15) 



1 

sin 9a(b) cos 9a(b) 
To diagonalize the spin Hamiltonian Eq. 2 in terms of the rotated spins S we introduce spin 
raising and lowering operators a+ and a on sublattice A by 



c+ _ 
^i£A - 


- S^^A + ^'S'jgA 


^i^A = 


- ^ieA ~ ''■^ieA 


^i&A = 


1 . + . 

= 2 - <«i 



(16) 

which obey the usual bosonic commutation relations in the large S limit. Similarly, operators 
b^ and b are introduced on sublattice B. After Fourier transformation this leads, up to a 
constant energy shift, to the linear spin wave Hamiltonian 



H 



SW — 2^ '[-^11 (o^k^k + a_k^-k 



H.3 bfbk + 6+ J-k ) + H. 



k"-k 



hi lOka+k + Oka-k 



HsA (b^b\ + 6k&-k) + ^31 (a^h + a\b^]^ + Ok^k + «-k&-k 



+^41 iatb\ + a:^k^k + «k&-k + a-k&: 



(17) 



neglecting higher order terms in a and b. Due to the two-sublattice structure the k-sum is 
restricted to half of the Brillouin zone, i.e. to momenta with cos(A;2:) + cos{ky) > 0. In the 
superfiuid and the supersolid phase the k-dependent coefficients in Eq. 17 are given by: 



9 



nil — Z— — T. H -n2i , -"33 — ^^ — 7, H -034 

sm 6 A sm 63 

i721 = ^^2Sin2^A7f , i/34=^^2sin2eB7['^ 
^3i = 7k^ -l-cos0ACOs^B + -Vism0^sin^B , H41 = H31 + 2-fl^^ (18) 

where 7j(. = |(cos(/i:2;) + cos{ky)) and 7,^ = cos(/c^) cos(/cy). The coefficients of the ffist 
order terms are required to vanish ||22| which leads to the conditions 

-— sin 6a = 2 cos 6a sin 6b + Vi cos 6b sin 6a + V2 cos 6*^ sin 6a 

TJ 

— ^ sin 6b = 2 cos 6'^ sin 6*^ + Vi cos 6'^ sin 6b + V2 cos 6*^ sin 6b ■ (19) 

These two equations determine the angles 6a and 6b for a given value of the magnetic field 
H^. The solutions of Eqs. 19 determine the phase diagram of the model. These equations 
fully coincide with the ones obtained by minimizing the free energy in the previous section. 
We have already used Eq. 19 to eliminate the magnetic field in the expressions for the 
coefficients Hn and if 33 in Eq. 18 of the spin wave Hamiltonian. However, for the Neel solid 
and the Mott insulator phase where both sin^^ = and sin^^ = the elimination is not 
possible and instead Hn and Hzz are given by 

Ell = Vi-V2- \h, H,, = Vi-V2 + \h, Neel Sohd 

Hii = Hs3 = -Vi -V2 + -H^ Mott Insulator (20) 

Eq. 17 is diagonalized by a generalized Bogoliubov transformation using the equation 
of motion idta\^ = [ak,Hsw]- with a^ oc e^"^''*. In the boson language the Bogoliubov 
transformation involves coupled density and phase modes. As a result we obtain the spin 
wave dispersion in the form 



^± 



(k) = Ij-f/'ii - if 21 + 2if|i - 2ii4i + iffg - if|4 ± [ i^Hf^ - H21 - i^la + H^^j + 

+4( [iin - ii2l] [iisi + H,i] + [ii33 + ii34] [ii31 - H21] ) ■ 

■( [ii33 - i^34] [ii31 + ^41] + [Hii + ii2l] [ii31 - ii2l] ) ) '^'} (21) 



10 



Typical dispersions are shown in Fig. 6 in the different phases and on the phase bound- 
aries. We now overview the dispersion relations in the four phases. 



i) In the Neel solid which is realized for V2 < Vi — 2 and Hz < 2J(Vi — V-i)^ — 4 the spin 
wave dispersion is given by 

o;±(k) = ^(^1 - V^Y - (27^^))' ± \Hz . (22) 

Thus, there are two excitation branches in a halved magnetic Brillouin zone. Both branches 
are gapped. 

a) In the superfluid there is a Goldstone mode of linear k dependence at small /c, and a 
well developed minimum around k* = (tt, tt). Taking the continuum limit carefully identifies 
this with the roton part of the helium dispersion. Explicitly, with s = sin^^ = sin^^, the 
dispersion in the extended zone is given by 

c.2(k) = 2 (1 - 7t^^) (2(1 - 7t'^) + s' [\/27f + (2 + V^W^^]) . (23) 

in) In the Mott insulator phase for fields Hz > 2(2 + Vi + V2) all spins are aligned along the 
magnetic field direction. There is a single gapped mode in the extended P* Brillouin zone 
with the dispersion given by 

.;(k) = ^-V,-V2- 27^'^ . (24) 

iv) Finally, in the supersolid phase one has a gapless linear mode, and a gapped one, again 
in the halved magnetic zone. 

To clarify the physics of the transitions we concentrate on the dispersion at k ^ and 
k ^ (tt, tt) at the phase boundaries. At the supersolid-Neel solid transition the critical mode 



is the Goldstone mode at small k. At the critical magnetic field, H^ = 2J(Vi — V'lf' — 4, 
which determines the Neel solid to supersolid boundary by the vanishing of the gap of the 
lower excitation branch of the solid, we perform the small k expansion for C(j_(k) from Eq. 22. 
For Vi > V2 + 2 where the solid exists at half-filling we obtain 

uj_ (k) ^ - ^ = k^ . (25) 

11 



This means that the hnear mode of the supersohd softens into a quadratic one at the 
boundary - signaUing the destruction of superfluidity - before hfting off into a gapped mode 
inside the sohd phase. This yields a quantum critical exponent z = 2. This value of z agrees 
with that of Chester [|ri|] and Cheng ||T^, but differs from that of Liu and Fisher [^], who 



obtain z = \. We feel, however, that the softening of the Goldstone mode is a physically 
reahstic picture, supporting our result. 

At the generic superfluid-to-Neel supersolid transition the critical mode is at k* = 
(7r,7r). Inside the superfiuid phase the roton minimum is at this wavevector. However in 
the solid, because of the zone-halving, this roton minimum is folded back to k = 0. In the 
superfiuid where sin 6'^ = sin 6b = s we study the small k expansion of the single mode 
in the neighborhood of k = (0,0) and k = (7r,7r). (Note that the Neel supersolid is only 
reahsed for V2 < l^i — 2.) 

uo^i-K, tt) - k) ^ 8A2 + (2 - 3A2 - V2S^) (26) 

where A^ = 2 + (s^/2) [V2 — 2 — V^i]. At the boundary to the Neel supersolid which is reached 



at a magnetic field H, = 2^(Vi - ¥2)^ - 4:{Vi + V2 + 2)/{Vi-V2+2) the mean field conditions 
in Eq. 19 tell that exactly at the transition the roton gap A disappears: the solidification is 
signalled by the softening-out of the roton mode of the superfiuid. The dispersion relation 
of the rotons also changes from a quadratic to a linear minimum, hence z = 1. 

Two remarks are in order here. First, recalling the original Landau argument about 
superfluidity it is clear that a vanishing roton energy leads to a vanishing critical velocity. 
So, while the superfluid order parameter remains finite through the supersolid transition, 
the critical velocity collapses to zero. Inside the supersolid phase it again assumes a finite 
value, as the second excitation branch becomes gapped. 

Second, one can raise the question of how this picture is going to be modified in the 
absence of an underlying lattice. In this continuum limit the modes which go soft are 
located at a finite magnitude of k, i.e. on a ring in momentum space. This means that 

12 



the phase space for these excitations is much larger than for the usual Goldstone modes, 
which are centered around A; ~ 0. It then is possible that these excitations may give rise to 
a fluctuation-induced first order transition instead of the second order one taking place on 



the lattice |^| 



Similar expansions can be used to study the case of half-filling. In this particle-hole 
symmetric case, not surprisingly both transitions have z = 1. In a recent Monte Carlo study 
IP the same z value was used in choosing the lattice size to study both the superfiuid- 
supersolid and supersolid-solid transitions, whereas we find z = 1 and z = 2, respectively 
off half-filling. It appears that our results call for the repetition of the numerical simulations 
with different z factors when p ^ 1/2. 

Finally at high fields, at the superfluid-to-Mott insulator transition the Goldstone mode 
softens out again, leading to z = 2, in agreement with earlier field theoretical predictions |]l| 
and numerical simulations p^ . 



We have repeated the spin wave calculation for the collinear ordering for an ordering 
wave vector k^, = (0, vr). In this case the coefficients of the Hamiltonian Eq. 17 outside the 
collinear solid and the Mott insulating phases are given by 



Hn={l + ^-^^] + 7 (Vi sin2 Or, - 2 - 2 cos^ 6^) cos(fc, 
^33 = ( 1 + ^^^] + \ {Vi sin^ 9r2 - 2 - 2 cos^ ^^2) cos(fc. 



sin 6* 



H21 = - sin^ Ori (2 + Vi) cos(A:^) , Hu = - sin^ ^^2 (2 + Vi) cos(A;^) , 

1 Vo (2) 

H^i = - [-2 - 2 cos 6ri cos 9r2 + Vi sin Ori sin dR2] cos{ky) + — sin Brx sin Br^ 7k 

i74i = Hzx + cos(A;j,) . (27) 

The MF conditions are read off from the vanishing of the terms linear in the spin wave 
operators as before 

Hz sin Qr\ = 2 cos 9Ri(sm 6ri + sin 6R2) + Vi sin 6ri (cos 6ri + cos OR2) + 2V2 cos 6r2 sin 6rx 

Hz sin 6r2 = 2 cos 9R2{sm 6ri + sin Or2) + Vi sin 6r2{co?, 6ri + cos 6*^2) + 2V2 cos Ori sin Or2 . (28) 

13 



For the superfluid and the Mott insulating phase the spin wave dispersions are obtained 
identical to the ones derived above in Eq. 23 and Eq. 24. In the coUinear solid with sin6'i^i = 
sin 6112 = the coefficients if n and H33 are replaced by 

Hn = V2- ^H, - cos(A;,) H^^ = V2 + ^H, - cos(A;,) . (29) 



The two gapped modes in the collinear solid for magnetic fields H^ < 2J(V2 — 1)^ — 1 
and V2 > 2 follow as 



a;±(k) = V(^2 - cos{K)Y - cos^iky) ± -if, . (30) 

As for the Neel supersolid, the collinear supersolid has one gapless linear mode at small 
k and a gapped one in the halved magnetic Brillouin zone which in the case of collinear 
ordering with wave vector k* = (0,7r) is determined by \ky\ < 7r/2. The transition from the 
superfiuid to the collinear supersolid is now driven by the softening of the roton mode at 
k* = (0, vr). The dynamical exponent is again z = 1. Also the exponents at the superfluid to 
collinear solid and at the solid to supersolid transition are identical to the exponents found 
for the Neel ordering transitions. 

IV. SIMULATIONS OF THE SOFT CORE MODEL 

Results at half— filling 

In this section we describe the results of numerical simulations, and compare them with 
the picture gained from the analytical considerations. Our Monte Carlo calculations are 
performed using a path integral representation on the BH partition function by discretizing 
the inverse temperature (3 into L^- intervals, (3 = Lt-At. A description of the technical 



details is contained in [0]. In order to characterize the phase diagram, we measure the boson 
winding number to determine the superfluid density ps- We also measure the density-density 
correlations c(l) and their Fourier transform, the structure factor S(k). 

c(I) = (n(j,r)n(j + l,r)) 
14 



^(k) = ^Ee'''(ca))- (31) 

jl 

Our normalization of the structure factor is such that if c(l) exhibits long range order, S'(k=f) 
will be proportional to the lattice volume A^ = L^, where L^ is the linear extent in the spatial 
dimension. If c(I) exhibits only short range order, S'(k^,) will be lattice size independent. 
Here k* = (tt, tt), (0, tt), (tt, 0) are the possible ordering wavevectors of the solid phase. 

At weak coupling or high temperatures, c(l) exhibits only short range order. For 1 small, 
c(l) is enhanced but very rapidly decays to its uncorrelated value p^. However, at low 
temperatures for sufficiently large interactions, the density-density correlations show long 
range oscillations. The associated structure factor 5'(k) evolves from being rather featureless 
to exhibiting a sharp peak at k* = (tt, tt) as Vi increases, and a peak at k* = (0, tt) or (vr, 0) 
as V2 increases. For our 2D system, for sufficiently large Vi, we expect a transition in the 
Ising universality class. That is, T,. is finite. In fact, if t = we have T^. = 0.567 Vi. But 
even for a zero temperature phase transition such as would occur at the Heisenberg point of 
the hard core model, one will still observe "long range order" at finite T when the diverging 
correlation length exceeds the spatial lattice size as T is lowered. In such instances, of course, 
a careful study of finite size effects is required to draw conclusions concerning the existence 
of long range order. Here we always report results for temperatures such that ^ > L^so that 
observables have taken on their ground state values. We have checked the scaling behavior 
to be sure that the ground state is genuinely ordered, when so claimed. 

Fig. 7 shows the superfiuid density Ps and structure factor S{7i, tt) as a function of Vi for 
Vo = 7 and V2 = 0. We see that at Vi ~ 2.5 there is a phase transition from a superfiuid to 
a solid phase. The transition on the 8x8 lattice shown is already rather sharp; finite size 
rounding in the raw data for the structure factor and superfiuid density near the transition 
point is further reduced as one goes to 10 x 10 lattices. That one has true diagonal long range 
order in the solid phase is confirmed by the fact that the structure factor scales linearly with 
the lattice volume. Indeed, at Vi = 8, S{7t, tt) is almost precisely 100/64 times as large on 
the 10 X 10 lattice than the 8x8. There does not appear to be any window of coexistence 
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between the superfluid and solid phases at half-filhng. To within hmits set by rounding, 
the transition points for S and p^ coincide almost precisely. We can make this statement 
more quantitative by performing the appropriate scaling analysis on the data. For example, 
we have plotted L'^S{n, vr) and L'^Ps versus Vi for different values of the exponent ratios 
a, b. Curves for different lattice sizes should cross at the same critical value of Vi for the 
appropriate choices of a, b. A complication is that the imaginary time lattice size must be 
scaled as the appropriate power of the spatial extent, and the dynamic exponent z could 
be different for the two transitions. Making the simplest assumption that z is the same, 
however, as was already suggested by the raw data, this scaling procedure shows that the 
transition points for the two observables are within 0.5% of each other. While the structure 
factors do indeed cross nicely, the superfluid density curves come together rather than pass 
through each other. This seems to be a rather generic feature of simulations of the Bose 
Hubbard model p5|] as opposed to related conserved current models [p^ ,p|. 



As V2 is increased, c(I) shows a similar transition from featureless uncorrelated behavior 
to long range order, although in this case V2 favors the formation of a "striped" collinear 
phase with alternating lines of occupied and empty sites. The structure factor S'(k) develops 
a peak at k,, = (vr, 0) or (0, vr). 

In order to determine whether V2 can drive a supersolid phase at half-filling, we turn 
on V2 close to the point where the transition between superfluid and solid occurs in Fig. 7. 
The density p = 0.5. We show in Fig. 8 a plot of ps and 5'(k) for k = (0,vr), (vr, 0) and 
(vr,vr). We see that V2 drives the Neel solid into a superfluid, and then at yet larger values 
causes the formation of a striped solid phase. Again, the plots suggest that there is no 
supersolid phase at p = 0.5. Scaling plots similar to those constructed at V2 = do not 
reveal any evidence for distinct critical points for superfluid and solid transitions to within 
our numerical accuracy. 

We can put data from Figs. 7 and 8 together with similar runs for different sweeps of Vi 
and V2 to obtain the ground state phase diagram of the soft core BH model at Vq = 7 and 
p = 0.5. This is shown in Fig. 9. At weak couplings we have a superfluid phase, while at 
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strong couplings there are two possible solids: checkerboard and striped. A strong coupling 
analysis predicts a phase boundary between the solid phases at V2 = \Vi The superfluid 
phase extends out along this line in a very robust manner, as opposed to the situation in 



ID, where the superfluid window was rather narrow. |^T| This is a consequence of the highly 
degenerate nature of the strong coupling (t = 0) ground state along the line Vi = 2V2. 
As can easily be seen, not only do the Neel and checkerboard solids have the same energy, 
but an infinite number of defect states are degenerate as well for Vi = V2. For example in 
a horizontally aligned collinear solid a whole column can be shifted up and down without 



energy cost. |l27|j28|l This large degeneracy stabilizes superfluidity, even at large coupling. 
We will comment further on this point when discussing the hard-core phase diagram. 
Results off half-filling 

Although it does not appear that the BH model exhibits a supersolid phase at p = 0.5, 
we can see the coexistence of diagonal and off-diagonal long range order when the filling 
is shifted away from p = 0.5. In Fig. 10 we show ps and ^(Tr, vr) for the same parameters 
as Fig. 7 except now p = 0.53. We see that although ps declines significantly when the 
solid forms, the excess boson density 6 = p — 0.5 (the magnetization m is spin language) 
remains mobile in the solid background. Indeed, simulations at different densities (we found 
supersolids out to dopings of 0.675) show that the tail in ps is precisely proportional to 6. 
Fig. 11 shows the analogous plot for a striped supersolid. Note that we have here separately 
displayed psx and psy. As expected, the superfluid density in the x and y direction is 
correlated with the direction in which the striped solid channels run, as determined by the 
ordering wavevector k^. = (tt, 0) or (0,7r). If we had separately measured p^a and pst on 
the two sublattices of the checkerboard solid, we would have found an analogous symmetry 
breaking. The nonzero value of psa — Psb is closely related to the appearance of a nonzero 
order parameter rrixa — ^^xb in the language of the spin Hamiltonian Eq. 2. 

If we were to use finite size scaling techniques to locate the precise phase boundaries, it 
would be necessary to scale the imaginary time length Lj. as a power of the spatial length 
L^, where z is the dynamic critical exponent. As we have earlier described, it may be that 
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different values of z are associated with the two transitions off half-filling, in which case the 
finite size scaling analysis is much more delicate. P,p!0| We do not see the necessity of such 
a study here, since the supersolid phase occupies an extended portion of the phase diagram, 
and its existence is not predicated on proving the distinctness of two transition points. 

Figures 10 and 11 provide compelling evidence for the existence of a supersolid phase. 
Our physical picture of this supersolid is one in which p = 0.5 of the bosons freeze into a 
rigid solid structure, while the remaining 5 remain mobile. As we have seen, a signal of long 
range order then is present in both the diagonal (niUj) and off-diagonal (aialj) channels. 

We have conducted our simulations of the BH Hamiltonian in the canonical ensemble, 
and have presented our results by specifying the density p rather than the chemical potential 
fi. In describing the nature of the phase diagram it is important to note that due to the 
existence of a gap in the solid phases, the /i-p relation is non-trivial. If the gap is nonzero, 
when we dope our system even slightly away from half-filling, the chemical potential is 
shifted by a considerable amount. In the language of the spin Hamiltonian, Eq. 2, a sizeable 
field Hz is required to change the magnetization of the gapped Ising phase. In Fig. 12 we 
illustrate this point by drawing the fi/Vi-1/Vi phase diagram. A sweep at constant chemical 
potential reveals a supersolid window. A sweep at fixed density skirts the pure solid and 
remains in the supersolid phase. This is why we see in Figs. 10 and 11 a supersolid for an 
extended region V > Krit rather than in some narrow region between phases exhibiting a 
single type of order. 

If we now examine densities p < 0.5, we find qualitatively similar results: a superfluid 
phase gives way to a striped supersolid phase as Vi increases. By these measures, hole 
or particle doping appears qualitatively similar. The same is true of the checkerboard 
supersolid, where results for hole doping are entirely reminiscent of the analogous particle 
doped case. 

In fact, however, something rather different does go on with particle and hole doping. 
In Fig. 13 we show the ground state energy as a function of doping for Vq = 7, Vi = 3, and 
V2 = 3. For these parameters, as we have seen, we have a striped supersolid off half-filling 



and a striped solid at p = 1/2. The change in slope of Eq at p = 1/2 reflects a jump 



in the chemical potential which is, in fact, just the gap in the solid phase. ^^. There is 
nothing particularly unusual here. The strange feature occurs for the checkerboard case. 
In Fig. 14 we show the ground state energy as a function of doping for Vq = 7, Vi = 3, 
and V2 = 0. The fact that Eq{p) is concave down for p < 0.5 indicates an instability to 
phase separation. Previous studies [^ have suggested the possibility of phase separation in 



systems with attractive boson interactions. However, we do not have these Lennard- Jones 
type potentials here, only purely repulsive ones. It is not immediately apparent why mobile 
holes (or particles) in a rigid solid background should segregate themselves. 

A possible explanation, however, is as follows: Consider an isolated doped hole in a 
checkerboard solid. In order to move to another site of the same sublattice, it must pass 
through an intermediate site on the opposite sublattice, a state of energy 2Vi. Thus the hole's 
effective hybridization is tefr = t'^/2Vi. (This sort of argument has previously been used to 
predict the shape of the phase boundary in the one dimensional extended BH Hamiltonian, 
in good agreement with simulations. ^T|) If two holes are near each other, the intermediate 
state is lower in energy, so the effective hybridization is increased. This suggests a possible 
mechanism for phase separation: increased mobility of holes which propagate coherently. Of 
course, the increase in tes is partially offset by the entropy cost of confining one hole near the 
other. Unfortunately, there appears to be an analogous increase in tes for doped particles 
which are proximate, so this reasoning does not explain the fact that E{p) is concave down 
for p < 0.5 only. Nevertheless, the simulations provide compelling evidence for a lack of 
particle-hole symmetry. 

In principle, one can also examine the issue of phase separation through anomalies in S'(k) 
for small k. However, our use of the canonical ensemble makes this approach non-trivial. 
Further work on the question of phase separation is needed. 
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V. SIMULATIONS OF THE HARD CORE MODEL 

We now examine the phase diagram in the hard core case. This is important to do 
for a number of reasons. First, it allows us to make a connection to the spin model limit, 
Eq. 2. Second, as we have seen at Vq = 7, some of the interesting transitions occur at Vi 
and V2 values which are getting rather large, while we expect in most physical situations 
that the on-site Vq should be substantially greater than the near neighbor interactions. One 
consequence of this, is that the doped bosons in the supersolid phase for our soft-core model 
could move on the occupied sublattice, since the cost of Vq was less than the coordination 
number z times the near neighbor interaction strengths. In the hard core model such multiple 
occupancies are forbidden, and we want to make sure that our conclusions are not affected 
by this change. 

Fig. 15 shows results for the superfluid density and structure factors for the half-filled 
case. We sweep V2 at fixed V^i = 3. A Neel phase appears at small V2. For larger V2 the 
superfluid phase appears before making a transition into a coUinear solid for yet larger V2. 
If Vi is sufficiently small, the Neel phase at weak V2 is eliminated, and the system remains 
superfluid down to V2 = 0. Data for this and other sweeps is summarized in Fig. 16 where 
the resulting ground state phase diagram is shown. Note that we find the superfluid-Neel 
solid transition at V2 = occurs at a value Vi close to 2t, which is the result expected based 
on the mapping to the spin model, Eq. 2. 

As in the soft core case the weak coupling superfiuid extends out along the V2 = V^i/2 
strong coupling boundary between the two solid phases. Unlike analogous studies [^ in 
ID, this superfiuid wedge is difficult to close, a phenomenon which we earlier explained by 
the large degeneracy of competing solid phases along the strong coupling line. We have 
conducted simulations along the line V2 = Vi /2 and find that the superfiuid density vanishes 
at Vi ~ 7. Interestingly, there is no inset of solid order at this point. This needs further 
study, for example to understand if some disordered dimer phase might exist in this regime, 
in analogy with related spin systems. 
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Fig. 17 is a plot for a doped lattice with S = p — 1/2 = 0.0625. The main difference is 
that, as in the soft-core case, there is a superfluid tail after the structure factor exhibits the 
transition into the solid phase. That is, there is a supersolid in the hard-core case as well. 
As expected, doping inhibits somewhat the formation of crystalline order, so that stronger 
couplings are required to induce the crystalline order as is seen by comparing the doped 
phase diagram. Fig. 18, with Fig. 16, the phase diagram at half-filling. Despite considerable 
rounding of the transitions, scaling analyses conclude that the regions where S is large are 
indeed ordered. 

Finally, Figs. 19a,b show the ground state energy as a function of filling for the hard-core 
model at Vi = 8, V2 = (Neel sohd), and V^i = 8, V2 = 4 (coUinear solid), respectively. The 
data are qualitatively similar to the soft-core model. In the collinear solid Eq is concave up, 
with a change in slope at p = 1/2 which is the gap. In the Neel solid Eq shows a tendency 
for phase separation. 

VI. RELATED ISSUES 

Up to now we have focussed on the ground state phase diagram of the BH Hamiltonian. 
It is interesting to consider also the behavior of the system at finite temperatures. Here the 
motion of doped bosons in the BH model which we have studied with our simulations has a 



close connection with the idea of "defectons" in a solid |jTj] where quantum tunneling caused 
by the finiteness of the de Boer parameter delocalizes lattice defects at low temperature. 
It is also of interest to study the behavior of the diffusion constant D for the full range of 
temperature. Here we expect that defects are localized at high T, and D first decreases 
exponentially as T is lowered in this classical regime. D should then exhibit a plateau as 
quantum diffusion takes over, and ultimately increase again as delocalization occurs. While 
they focus largely on the behavior of single defectons, Andreev and Lifshitz also consider 
the possibility of long range Coulomb interactions causing localization into a "defecton 
superlattice" . Our insulating checkerboard solid is in fact an illustration of this. The Bose- 
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Hubbard model with only on-site Vq, has no solid phase at p = 0.5, but when Vi is turned 
on, an ordered lattice does form. 

We have focussed here on zero temperature, the finite temperature phase diagram of the 
2D BH model would be interesting to study as well. The solid transitions are in the Ising 
universality class, and hence have a finite T^. Similarly one expects a Kosterlitz-Thouless 
type finite critical temperature for the superfiuid transition. As for the topology of the phase 
diagram, several possibilities have been explored by Liu and Fisher P| . One intriguing case is 
the appearence of a tetracritical point; where the three ordered phases (superfiuid, supersolid 
and solid) come together, giving way to the disordered phase with further increase of the 
temperature. This happens within a limited, but finite range of parameters on the mean 



field level. The corresponding scaling theory was developed by Nelson and Fisher ||33|. Other 
alternatives include a supersolid phase which exists only at finite temperatures, and that 
the tetracritical point is split into bicritical points. We hope to take up some of the issues 
in a further publication. 

In the path integral representation of the BH partition function used in our simulations, 
the particle number conservation leads to boson "world lines" propagating in the original 2 
spatial dimensions plus an additional imaginary time direction which runs from to (3. This 
picture has been used to suggest close analogies between the physics of vortices in Type II 
superconductors [^ and the phase diagram of the 2D Bose-Hubbard Hamiltonian. Frey, 
Nelson, and Fisher, [^ have recently discussed both thermally driven and quantum phase 



transitions, for example as caused by the introduction of defects or interstitials into the 
Abrikosov lattice, in these vortex systems. This also has close connections with the results 
we have discussed here. 

VII. CONCLUSIONS 

In this paper we have considered quantum phase transitions in the Bose-Hubbard hamil- 
tonian. We identified several phases: solid and supersolid phases with Neel and coUinear 
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patterns, furthermore a superfluid and a Mott type insulating phase. The phase diagram 
has been determined analytically and the spin-wave spectrum has been calculated. The dy- 
namical critical exponents at each transitions were calculated and preexisting controversies 
were settled. Our numerical work - utilizing Quantum Monte Carlo methods - provided 
a detailed study of the different phases. Concerning the phase diagram the existence of 
supersolid phases has been forcefully confirmed. These phases exist only off half-filling, in 
accordance with the mean field results, but in disagreement with some recent claims. The 
possibility of phase separation in the model has been investigated as well, and provides 
evidence for a violation of the previously assumed particle-hole symmetry of the model. 
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Figure Captions 

Fig. 1: Mean field phases MP of tlie XXZ spin Hamiltonian Eq. 2 on a 2D square lattice. 

Fig. 2: Hard core mean field phase diagram at half-filling p = | from comparing the energies 

of superfluid, Neel and collinear solid. 

Fig. 3: Hard core mean field phase diagram (bold lines) away from half-filling for m = 

|p — || = 0.2 from comparing the energies of superfluid, Neel and collinear supersolid. Thin 

lines indicate the phase boundaries at half-filling m = (see Fig. 2). 

Fig. 4: Hard core mean field phase diagram, magnetisation m versus V2, for fixed Vi = 3. 

SS denotes the supersolid phases. 

Fig. 5: Hard core mean field phase diagram, magnetisation m versus Vi, for fixed V2 = 1. 

Fig. 6: Spin wave dispersions in the (a) Neel solid, (b) at the Neel solid-Neel supersolid 

boundary, (c) in the Neel supersolid, (d) at the Neel supersolid-superfluid boundary, and 

(e) in the superfluid. In all plots V2 and Vi are flxed to V2 = 1.5, Vi = 4. and the magnetic 

fleld h {Hz in the text) is varied. 

Fig. 7: The superfluid density ps and S'(7r, tt) as a function of Vi for Vq = 7 and V2 = 0. 

The density p = 0.5 and (3 = 4. The transitions in ps and 5'(7r, tt) appear to occur at roughly 

the same value of Vi. 

Fig. 8: The superfluid density and structure factor as a function of V2 at Vq = 7, Vi = 2.75, 

and p = 0.5. 

Fig. 9: The ground state phase diagram of the BH model at p = 0.5 and with a soft core 

on-site repulsion Vq = 7. 

Fig. 10: The superfluid density and structure factor for the same parameters as in Fig. 7, 

except now the system is doped to p = 0.53. A superfluid tail remains in the (checkerboard) 

solid phase. 

Fig. 11: The superfluid density and structure factor for the same parameters as in Fig. 8, 

except now the system is doped to p = 0.56. A superfluid tail remains in the (striped) solid 

phase. 
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Fig. 12: The qualitative T = pliase diagram of the Bose-Hubbard model is illustrated. 
A sweep at constant /i (full arrow) could cut across the phase boundaries as shown, reveal- 
ing a supersolid window, while a sweep of constant density (dashed arrow) remains in the 
supersolid phase at strong coupling. The Mott insulating phase at large fj, has density one 
boson per site. 

Fig. 13: The ground state energy as a function of density for Vq = 7, Vi = 3, V2 = 3. 
Fig. 14: The ground state energy as a function of density for Vq = 7, Vl = 3, V2 = 0. 
Fig. 15: Superfiuid density and structure factors versus V2 for Vi = 3. The density p = 
0.500. 

Fig. 16: The phase diagram of the half-filled hard-core model. The dashed lines are the 
results of the Mean Field analysis presented earlier in the paper. 

Fig. 17: Superfiuid density and structure factors versus V2 for Vi = 3. The density p = 
0.563. 

Fig. 18: The phase diagram of the doped hard-core model. 

Fig. 19: Ground state energy versus density for the hardcore model, (a) Vi = 4, V2 = 
and (b) Vi = 4, 1/2 = 4. 
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